% This file generates the impulse response presented in Figure C1 in 
% Appendix C of the paper. It is being called by the file 
% 'FIGs_C1_C2_C3_C4_Medium_Scale_Main.m'. This file runs the Dynare file 
% 'FIG_C1_base_model.mod'.

close all; clear all; clc;

global y0 long drop replic iorder prune p_int peg sug dr ga alpha gi mcs ...
       rks xiw sigma pistar gammaw delta eta b beta ys chi phi ylsf klsf ...
       eta beta b govs delta gup gi sigma wrsf wsf chi xiwf pistar gammaw prune

%% Options:

fig_plot        =   0;
prune           =   1;          % set to 1 for pruning
third           =   0;          % set to 1 to do a third order approximation, 0 to do second
growth_source   =   1;          % set to 0 for both A and I, 1 for just A, 2 for just I

replic          =   50;         % number of replications when calculating irfs
drop            =   0;          % number of periods to drop when calculating irds
T               =   500;        % number of periods to simulate
burn            =   100;        % number of periods to burn-in after steady state
peg             =   8;          % number of periods to peg
long            =   22;         % length of IRFs

alpha           =   1/3;
phi             =   0.001;
gad             =   1.0033;
gid             =   1.000;
thetag          =   0.5;        % parameter governing utility from government spending
gad             =   1;

iorder          =   1;

e11              =   clock;

%load estimated_params

psi1            =   1/((1-phi)*(1-alpha));
psi2            =   (alpha*(1-phi))/((1-phi)*(1-alpha));
check1          =   psi1*(phi-1 + alpha*(1-phi)) + 1;
check2          =   psi2*(phi-1 + alpha*(1-phi)) + alpha*(1-phi);

omegag          =   0.0;
beta            =   0.995;
theta           =   11;
sigma           =   11;
pistar          =   1.005;
pistar          =   1;
phi_fix         =   0.00;
% Desired output growth rate:
gupd            =   gad^(1/((1-phi_fix)*(1-alpha)))*gid^(alpha/(1-alpha));
gamma2          =   0.05;
delta           =   0.025;
xiw             =   0.66;
xip             =   0.66;
b               =   0.7;
chi             =   1;
kappa           =   4;
rhoA            =   0.95;
gammap          =   0;
gammaw          =   0;
alphapi         =   1.5;
alphay          =   0.2;
rhoi            =   0.8;
rhoinvs         =   0.8;
rhob            =   0.6;
rhoL            =   0.8;
suA             =   0.01;
suL             =   0.01;
sug             =   0.01;
rhog            =   0.9;
sub             =   0.01;
sur             =   0.0025;
eta             =   6;
suinvs          =   0.04;
alphax          =   0;

if growth_source == 0
    ga              =   gad^(1-phi);
    gi              =   gid;
    gup             =   ga^(1/((1-phi)*(1-alpha)))*gi^(alpha/(1-alpha));
else
    if growth_source == 1
        ga          =   gupd^((1-phi)*(1-alpha));
        gi          =   1;
        gup         =   ga^(1/((1-phi)*(1-alpha)))*gi^(alpha/(1-alpha));
    else
        ga          =   1;
        gi          =   gupd^((1-alpha)/alpha);
        gup         =   ga^(1/((1-phi)*(1-alpha)))*gi^(alpha/(1-alpha));
    end
end

%% Compute steady states.

As              =   1;

ps              =   ((1-xip*pistar^((gammap-1)*(1-theta)))/(1-xip))^(1/(1-theta));
ss              =   (1-xip*pistar^(-theta*(gammap-1)))^(-1)*((1-xip)*ps^(-theta));
ints            =   pistar*gup/beta - 1;
mcs             =   ((theta-1)/theta)*((1-xip*pistar^((gammap-1)...
                    *(1-theta)))/(1-xip))^(1/(1-theta))*(1-xip*beta*pistar^((gammap-1)...
                    *(1-theta)))^(-1)*(1-xip*beta*pistar^(-theta*(gammap-1)));
rks             =   gup*gi/beta - (1-delta);
gamma1          =   rks;
ws              =   ((As*mcs*rks^(alpha*(phi-1)))/(phi^(-phi)*(1-phi)^(phi-1)...
                    *(alpha^(-alpha)*(1-alpha)^(alpha-1))^(1-phi)))^(1/((1-alpha)*(1-phi)));
wrs             =   ((1-xiw*gup^(sigma-1)*pistar^((gammaw-1)*(1-sigma)))/(1-xiw))^(1/(1-sigma))*ws;
kls             =   (alpha/(1-alpha))*gup*gi*(ws/rks);
gks             =   (phi/(1-phi))*(1/alpha)*gup^(-1)*gi^(-1)*rks;
xls             =   As*mcs*gks^(phi)*kls^(phi*(1-alpha) + alpha)...
                    *gup^(alpha*(phi-1))*gi^(alpha*(phi-1));
xls2            =   As*(1/ss)*gks^(phi)*kls^(phi*(1-alpha) + alpha)...
                    *gup^(alpha*(phi-1))*gi^(alpha*(phi-1));
yls             =   xls - gks*kls;
cls             =   (1-omegag)*yls - (1-(1-delta)*gup^(-1)*gi^(-1))*kls;
ls              =   (((sigma-1)/sigma)*(cls)^(-1)*(1/eta)*(((gup-b) - beta*b...
                    *(1-b*gup^(-1)))/((1-b*gup^(-1))*(gup-b)))*(wrs/ws)^(sigma*chi)...
                    *wrs*(1-beta*xiw*gup^(sigma-1)*pistar^((sigma-1)*...
                    (1-gammaw)))^(-1)*(1-beta*xiw*gup^(sigma*(1+chi))...
                    *pistar^(sigma*(1+chi)*(1-gammaw))))^(1/(1+chi));
ys              =   yls*ls;
cs              =   cls*ls;
ks              =   kls*ls;
gs              =   gks*ks;
xs              =   xls*ls;
govs            =   omegag*ys;
invss           =   (1-(1-delta)*gup^(-1)*gi^(-1))*ks;
F               =   ((1-mcs*ss)/mcs)*xs;
% SS MUC:
lams            =   (1/cs)*((gup - b) - beta*b*(1-b*gup^(-1)))/((1-b*gup^(-1))*(gup-b)); 
x1s             =   (lams*xs*mcs)/(1-xip*beta*pistar^(-theta*(gammap-1)));
x2s             =   (lams*xs)/(1-xip*beta*pistar^((gammap-1)*(1-theta)));
f1s             =   (1-beta*xiw*gup^(sigma*(1+chi))*pistar^(sigma*(1+chi)*(1-gammaw)))^(-1)...
                    *(eta*(ws/wrs)^(sigma*(1+chi))*ls^(1+chi));
f2s             =   (1-beta*xiw*gup^(sigma-1)*pistar^((sigma-1)*(1-gammaw)))^(-1)...
                    *(lams*(ws/wrs)^(sigma)*ls);
vws             =   (1-xiw*(gup^(-1)*pistar^(gammaw-1))^(-sigma*(1+chi)))^(-1)...
                    *((1-xiw)*(wrs/ws)^(-sigma*(1+chi)));
% Check to make sure capital works out correctly:
k_check         =   ks - gi*gup*(alpha*(1-phi)*(mcs/rks)*(ss*xs + F)); 
% Check to make sure labor works out correctly:
l_check         =   ls - (1-alpha)*(1-phi)*(mcs/ws)*(ss*xs + F); 
% Check to make sure gamma works out:
gamma_check     =   gs - phi*mcs*(ss*xs + F); 
profit_check    =   xs - ws*ls - rks*ks*gup^(-1)*gi^(-1) - gs;
reset_price_check = ps - (theta/(theta-1))*x1s/x2s;
reset_wage_check =  wrs - (sigma/(sigma - 1))*f1s/f2s;

%% Find flexible price equilibrium.
xipf            =   0.00;
xiwf            =   0.00;

% Compute steady states. 
psf             =   ((1-xipf*pistar^((gammap-1)*(1-theta)))/(1-xipf))^(1/(1-theta));
ssf             =   (1-xipf*pistar^(-theta*(gammap-1)))^(-1)*((1-xipf)*psf^(-theta));
intsf           =   pistar*gup/beta - 1;
mcsf            =   ((theta-1)/theta)*((1-xipf*pistar^((gammap-1)...
                    *(1-theta)))/(1-xipf))^(1/(1-theta))*(1-xipf*beta...
                    *pistar^((gammap-1)*(1-theta)))^(-1)*(1-xipf*beta*pistar^(-theta*(gammap-1)));
rksf            =   gup*gi/beta - (1-delta);
wsf             =   ((mcsf*rksf^(alpha*(phi-1)))/(phi^(-phi)*(1-phi)^(phi-1)...
                    *(alpha^(-alpha)*(1-alpha)^(alpha-1))^(1-phi)))^(1/((1-alpha)*(1-phi)));
wrsf            =   ((1-xiwf*gup^(sigma-1)*pistar^((gammaw-1)*(1-sigma)))/(1-xiwf))^(1/(1-sigma))*wsf;
klsf            =   (alpha/(1-alpha))*gup*gi*(wsf/rksf);
gksf            =   (phi/(1-phi))*(1/alpha)*gup^(-1)*gi^(-1)*rksf;
xlsf            =   mcsf*gksf^(phi)*klsf^(phi*(1-alpha) + alpha)...
                    *gup^(alpha*(phi-1))*gi^(alpha*(phi-1));
xls2f           = (1/ss)*gksf^(phi)*klsf^(phi*(1-alpha) + alpha)...
                    *gup^(alpha*(phi-1))*gi^(alpha*(phi-1));
ylsf            =   xlsf - gksf*klsf;

% Numerically solve for C/L, G/L, and L.
options         =   optimset('TolX',1000e-25,'TolFun',1000e-25,...
                    'MaxFunEvals',1000e+25,'MaxIter',50,'LargeScale','off','display','iter');
guess           =   0.33;                   % Guess lsf.
[ll,fval]       =   fsolve(@FIGs_C1_C2_C3_C4_find_flex_price,guess,options);
lsf             =   ll;
clsf            =   ylsf - govs/lsf - (1-(1-delta)*gup^(-1)*gi^(-1))*klsf;
lsf_check       =   (((sigma-1)/sigma)*(clsf)^(-1)*(1/eta)*(((gup-b)...
                    - beta*b*(1-b*gup^(-1)))/((1-b*gup^(-1))*(gup-b)))...
                    *(wrsf/wsf)^(sigma*chi)*wrsf*(1-beta*xiwf*gup^(sigma-1)...
                    *pistar^((sigma-1)*(1-gammaw)))^(-1)*(1-beta*xiwf*gup^(sigma*(1+chi))...
                    *pistar^(sigma*(1+chi)*(1-gammaw))))^(1/(1+chi));
ysf             =   ylsf*lsf;
csf             =   clsf*lsf;
ksf             =   klsf*lsf;
gsf             =   gksf*ksf;
xsf             =   xlsf*lsf;
govsf           =   omegag*ysf;
invssf          =   (1-(1-delta)*gup^(-1)*gi^(-1))*ksf;
Ff              =   ((1-mcsf*ssf)/mcsf)*xsf;
% SS MUC:
lamsf           =   (1/csf)*((gup - b) - beta*b*(1-b*gup^(-1)))/((1-b*gup^(-1))*(gup-b)); 
x1sf            =   (lamsf*xsf*mcsf)/(1-xipf*beta*pistar^(-theta*(gammap-1)));
x2sf            =   (lamsf*xsf)/(1-xipf*beta*pistar^((gammap-1)*(1-theta)));
f1sf            =   (1-beta*xiwf*gup^(sigma*(1+chi))*pistar^(sigma*(1+chi)...
                    *(1-gammaw)))^(-1)*(eta*(wsf/wrsf)^(sigma*(1+chi))*lsf^(1+chi));
f2sf            =   (1-beta*xiwf*gup^(sigma-1)*pistar^((sigma-1)*(1-gammaw)))^(-1)...
                    *(lamsf*(wsf/wrsf)^(sigma)*lsf);
% Check to make sure capital works out correctly:
k_checkf        =   ksf - gi*gup*(alpha*(1-phi)*(mcsf/rksf)*(ssf*xsf + Ff)); 
% Check to make sure labor works out correctly:
l_checkf        =   lsf - (1-alpha)*(1-phi)*(mcsf/wsf)*(ssf*xsf + Ff); 
% Check to make sure gamma works out:
gamma_checkf    = gsf - phi*mcsf*(ssf*xsf + Ff); 
reset_price_checkf = psf - (theta/(theta-1))*x1sf/x2sf;
reset_wage_checkf  = wrsf - (sigma/(sigma - 1))*f1sf/f2sf;
profit_checkf 	=   xsf - ws*lsf - rks*ksf*gup^(-1)*gi^(-1) - gsf;
vwsf            =   (1-xiwf*(gup^(-1)*pistar^(gammaw-1))^(-sigma*(1+chi)))^(-1)...
                    *((1-xiwf)*(wrsf/wsf)^(-sigma*(1+chi)));
gaps            =   ys/ysf;

%% Shocks, model solution, and variable positions. 

rhor            =   0.00;
css             =   cs;
sii             =   0.002;          % Standard deviation of anticipated shock.

save parameter_main_gov_peg beta b gamma1 gamma2 kappa delta pistar sigma...
     chi eta xiw gammaw alpha phi F theta xip gammap alphapi alphay rhoi...
     rhob rhoL rhoA rhor rhoinvs sub suL suA sur suinvs ints ys gup gi lams...
     ls ks ws wrs ps gs x1s x2s f1s f2s xs mcs ss As vws css lamsf csf wrsf ...
     f1sf f2sf wsf lsf mcsf ssf xsf gsf psf x1sf x2sf ysf ksf invssf xipf ...
     xiwf Ff gaps alphax omegag rhog sug thetag govs

dynare FIG_C1_base_model noclearall nolog nostrict

dr              =   oo_.dr;
y0              =   oo_.dr.ys;

% Variable positions:
p_lam           =   1;
p_c             =   2;
p_rk            =   3;
p_z             =   4;
p_mu            =   5;
p_invs          =   6;
p_int           =   7;
p_infl          =   8;
p_wr            =   9;
p_f1            =   10;
p_f2            =   11;
p_khat          =   12;
p_v             =   13;
p_s             =   14;   
p_L             =   15;
p_pr            =   16;
p_x1            =   17;
p_x2            =   18;
p_w             =   19;
p_y             =   20;
p_k             =   21;
p_eb            =   22;
p_eL            =   23;
p_A             =   24;
p_er            =   25;
p_einvs         =   26;
p_dy            =   27;
p_dc            =   28;
p_dI            =   29;
p_dL            =   30;
p_dw            =   31;

guess           =   zeros(1,peg);
[xx,fval]       =   fsolve(@FIGs_C1_C2_C3_C4_find_peg_shocks,guess);
e1              =   [0 0 0.01 xx(1) 0 0 xx(2) xx(3) xx(4) xx(5) xx(6) xx(7) xx(8)]';

%% Impulse Response Functions. 

% Modification of Dynare's IRF command:
irf_ss          =   FIGs_C1_C2_C3_C4_irf_dynare(dr, e1, long+1, drop, replic, iorder);
r_irf_ss        =   zeros(1,long);

for jx = 1:long
    r_irf_ss(jx)    =   irf_ss(p_int,jx) - irf_ss(p_infl,jx+1);
end

% Now do it without the ZLB:
xx              =   zeros(1,peg);
e2              =   [0 0 0.01 xx(1) 0 0 xx(2) xx(3) xx(4) xx(5) xx(6) xx(7) xx(8)]';
irf2            =   FIGs_C1_C2_C3_C4_irf_dynare(dr, e2, long+1, drop, replic, iorder);
r_irf2          =   zeros(1,long);

for jx = 1:long
    r_irf2(jx) = irf2(p_int,jx) - irf2(p_infl,jx+1);
end

e22             =   clock;

elapsed         =   etime(e22,e11)/60
